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We investigate tree tensor network states for quantum chemistry. Tree tensor network states 
represent one of the simplest generalizations of matrix product states and the density matrix renor- 
mahzation group. While matrix product states encode a one-dimensional entanglement structure, 
tree tensor network states encode a tree entanglement structure, allowing for a more flexible descrip- 
tion of general molecules. We describe an optimal tree tensor network state algorithm for quantum 
chemistry. We introduce the concept of half-renormalization which greatly improves the efficiency 
of the calculations. Using our efficient formulation we demonstrate the strengths and weaknesses 
of tree tensor network states versus matrix product states. We carry out benchmark calculations 
both on tree systems (hydrogen trees and vr-conjugated dendrimers) as well as non-tree molecules 
(hydrogen chains, nitrogen dimer, and chromium dimer). In general, tree tensor network states 
require much fewer renormalized states to achieve the same accuracy as matrix product states. In 
non-tree molecules, whether this translates into a computational savings is system dependent, due 
to the higher prefactor and computational scaling associated with tree algorithms. In tree like 
molecules, tree network states are easily superior to matrix product states. As an illustration, our 
largest dendrimer calculation with tree tensor network states correlates 110 electrons in 110 active 
orbit als. 



INTRODUCTION 

Currently, there is much effort devoted to finding effi- 
cient numerical techniques for strongly correlated elec- 
trons. Amongst several approaches, the Density Ma- 
trix Renormahzation Group (DMRG)[l|, Q has pro- 
vided many new insights in challenging systems. In re- 
cent years, efficient DMRG implementations have also 
appeared for quantum chemistry. [3|49| Originally, the 
DMRG was formulated as an algorithm in the language 
of Wilson's numerical renormahzation group. However, 
more recently, attention has expanded to the underlying 
class of wavefunctions optimized by the DMRG, which 
are the Matrix Product States (MPS). 

MPS provide a compact description of entanglement 
in ID systems, thus DMRG calculations for chain-like 
molecules are very efficient. [Tol, [ll] However, when used 
in 2D and 3D systems, much larger bond-dimensions, 
usually denoted by M and referred to as the number 
of renormalized states in DMRG calculations, must be 
used to reach a target accuracy. The need for larger 
M reflects the sub-optimal representation of 2D and 3D 
entanglement by the MPS wavefunction structure. 

The generalization of MPS to tensor network states 
(TNS) provides a natural way to compactly describe 2D 
and 3D entanglement. There are several families of TNS 
that differ in the way the entanglement is encoded, and a 
brief overview of the classes of TNS is given in Refs. [12I - 
171 . While the formal properties of general TNS are ap- 
peahng, efficient computation with these states lags far 
behind computation with MPS. In this work, we explore 
efficient computation for quantum chemistry with Tree 
Tensor Network States (TTNS),pJ, llsHH one of the 
simplest families of TNS. 



TTNS encode a tree entanglement structure, as illus- 
trated in Figure [H A tree of maximal degree Z has 
at most Z neighbours at any site. MPS are a special 
case of TTNS with Z = 2. The absence of cycles in 
a TTNS simplifies many algorithms. In particular, the 
variational minimization of TTNS wavefunctions can be 
carried out by generalizing the DMRG algorithm used 
for MPS, from a two-block formulation, to a Z-block 
formulation. Multi-block DMRG algorithms have been 
studied for some time422r.27j although to the best of our 
knowledge, their interpretation in terms of the under- 
lyin g cl ass of TTNS variational states first appeared in 
Ref. 119|. Thus, Ref. 19 can be considered as a starting- 
point for the current work. 

In Ref. 19 the authors considered a prototype applica- 
tion of TTNS to quantum chemistry with a minimal basis 
beryhium atom calculation. Although this provided evi- 
dence that for a given M, TTNS of degree Z > 2 capture 
more entanglement than the corresponding MPS with the 
same M, for practical computation many questions were 
left unanswered. For example, although TTNS are more 
flexible than MPS for describing entanglement, there is 
an implicit trade-off of higher computational scaling with 
M. Furthermore, the calculation in Ref. 19 used an un- 
realistically small M = 3, whereas typical DMRG calcu- 
lations use M = 1000 - 10000. Thus, the relevance of 
TTNS for quantum chemistry calculations remains to be 
established. 

In the current work, we address this question by de- 
scribing and implementing an efficient TTNS algorithm 
that can be used for realistic calculations. Our algorithm 
is formulated to have an optimal computational scaling 
for quantum chemistry Hamiltonians. We achieve a much 
lower cost than earlier formulations through a "half"- 
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FIG. 1. Examples of trees. The left panel shows a tree with 
degree Z = 3 and depth A = 4, and the right panel shows a 
tree with degree Z — \ and depth A = 3. 



renormalization transformation which exactly maps the 
multi-block DMRG to a conventional, efficient, two-block 
DMRG. We also address the issue of orbital ordering on 
trees which is necessary to use TTNS in chemistry. Using 
our efficient implementation, we assess the performance 
of TTNS relative to MPS used in quantum chemistry 
DMRG calculations. We compare TTNS and MPS on 
several benchmark molecular systems, including ideal hy- 
drogen trees and chains, and benchmark molecules such 
as the nitrogen and chromium dimers. Finally, to demon- 
strate the power of TTNS, we carry out complete active 
space calculations in 7r-conjugated dendrimer molecules, 
correlating up to 110 electrons in 110 orbitals. 



OVERVIEW OF THE DMRG ALGORITHM 
BASED ON MPS 

We first present a brief overview of MPS wavefunctions 
and the DMRG algorithm that optimizes their energy. 
This will allow us to establish basic notation which will 
be used to discuss TTNS in the next section. However, as 
the presentation is not entirely self-contained, for further 
details we refer to additional articles and reviews. jl|, 0, 
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We use the following notation for MPS: is the many- 
body basis at site i, k is number of sites, and M is 
the number of renormalized states. The (one-site) MPS 
wavefunction is obtained by expressing the coefficient of 
the determinant |ni...n/c) as a product of matrices for 
each occupancy ni...n/c, 

1^) = ^ L^\..L^^-^V'^^R^^+^..R^'^ |ni...n/e) (1) 



For an MPS with M renormalized states, the matrices 
are of maximum dimension M x except for the first 
and last, which are of maximum dimension 1 x M and 
M X 1 respectively. Note that the MPS is invariant to a 
number of transformations of the matrices, [sl. [29j We re- 
move this invar iance by choosing an MPS canonical form. 
In the canonical form at site i, the rotation matrices to 



the left of site i are constrained to satisfy orthonormality 
conditions L^^^'^'L"^* = 1, while those to the right of 
site i are constrained to satisfy R"^*!^"^*^ = 1- 

The left and right rotation matrices allow us to define 
left and right renormalized many-body states, and 
Ir^), respectively. These renormalized representations 
are used to construct the computational intermediates 
(renormalized operators) in the DMRG algorithm. Car- 
rying out the matrix multiplications from L^^ through 
L^*-i, and from R"^'^ through R^*+i, we obtain 



\k- 



E 

ni . ..rii- 



\ni...ni-i) 



J2 R"'+^..R"Mni+i-nfc). 



(2) 



(3) 



The orthonormality conditions on the rotation matrices 
L and R imply that the renormalized bases are orthonor- 
mal 



(4) 
(5) 



The DMRG (canonical) form of the MPS wavefunction 
is obtained by rewriting the MPS wavefunction ([T]) in 
terms of the renormalized many-body bases and 



I*) 



E 

li — 1 '^i l^i 



(6) 



In this interpretation, i/^^* is viewed as a wavefunction 
coefficient vector in a Hilbert space spanned by the renor- 
malized product states, \li-\niTi) = \ni) \ri). 

Computations using MPS involve tensor- tensor con- 
tractions. To express such operations, it is helpful to 
use a graphical representation, shown in Figure [21 Each 
vertex is a tensor and the number of edges connected to 
the vertex determines the tensor rank. In the case of an 
MPS, each L^* or R^* is a rank-3 tensor, represented by 
a vertex with three edges. Here, we always choose the 
vertical index to represent n^. The MPS wavefunction 
is obtained by contracting the horizontal edges of all the 
tensors (Fig.[2](B)), leading to Eq. ([T|). The computation 
of the overlap of two MPS is shown in Fig. [2] (C). 

We optimize the MPS energy, by minimizing the La- 
grangian (^|^|^) - (A(^|^) - 1) with respect to the 
tensors in the MPS. In the (one-site) DMRG sweep al- 
gorithm, this minimization is carried out with respect to 
a single tensor at a time. In step i of the DMRG sweep, 
the MPS is expressed in the DMRG form (|6]), and the 
coefficient vector i/? is optimized, holding the rotation 
matrices to the left and right of the site which define the 
bases Ir^), fixed. Because the energy is a quadratic 

form in '0, minimization leads to a standard eigenvalue 
problem, 

^ {k.^n,n\H\li_,n',r[) ^l^^^ - = (7) 
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FIG. 2. Graphical representations of MPS. (A) one-site coefficient vector '0^* is a rank-3 tensor, (B) the MPS wavefunction is 
given by contracting all horizontal bonds, (C) the overlap of two MPS wavefunctions can be efficiently computed recursively. 
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FIG. 3. Graphical representations of DMRG wavefunction. 
Canonical form of MPS wavefunction (top) can be re- written 
as a block diagram in DMRG language (bottom). 

In subsequent steps of the sweep, the MPS is transformed 
to the DMRG form at successive sites, and the coefficient 
vectors at these sites are optimized. 

The most expensive operation in the sweep is com- 
puting and performing Hip to solve the eigen- 
value problem ([7j). This is because the second-quantized 
Hamiltonian contains a large number of terms, 

H = ^ hija\aj + ^ X] Vijika\a}-akai' (8) 

ij ijkl 

and for each term we need its matrix representation in 
the basis {Ih-iriiri)}. There are two generic strate- 
gies to handle the large number of terms. The first, 
used in existing quantum chemistry DMRG irnplemen- 
tations, is through Complementary Operators.^^ 0, [30| 
The second uses the more recent concept of Matrix Prod- 
uct Operators.[31-35] 

Complementary operators are a way to maximize the 
reuse of intermediates. For example, there are 0{k^) 
terms in the summation (|8j), and each expectation value 
individually is of 0{M^) cost, leading to a naive scal- 
ing of 0{M^k'^) for the energy. However, much infor- 
mation can be reused between terms. For example, the 
two terms {ala^asa^) and (alalaQaj) involve the same 
partial expectation value over a|a2. Complementary op- 
erators reuse and combine such partial traces. At site 
i in the DMRG sweep, we partition the Hilbert space 
into two subspaces: A, containing the left block of sites 
(sites — 1) and site i, spanned by renormalized 



states \li-ini)^ and B, containing the right block of sites 
i + 1 . . . /c, and spanned by renormalized states \ri). (See 
Fig. [3j). H is correspondingly partitioned as: 

H = Ha^Hb^ Hab (9) 

Ha and Hb are act locally on A and B, respectively, and 
have non-trivial expectation values with only \li-ini) and 
\ri) separately. Hab describes the interactions between 
A and B, and is given by a sum of products acting sepa- 
rately on the two spaces 

Hab= E {afsr+afst + afam) 

^^A,eB ^^^^ 

+ ^ (af a/f^f t + af^afQf^) + adjoint. 

i>jeA 

where Pij , Qij , and Si are the Complementary Operators^ 

= ^ VijikCLkai (11) 
kieB 

Qfj = ^ {vikji - Vikij)a\ai (12) 

kleB 

Sf = ^tijCtj + ^ Vijika]dkai. (13) 
jeB jkieB 

Using complementary operators, the total complexity of 
evaluating and performing H"^ is reduced to 

0{M^k^ + M'^k'^) per sweep, including the cost of con- 
structing the components of H in the partitioned form 
(renormalization) . y This is the standard procedure to 
evaluate the energy in the DMRG algorithm. 

An alternative approach to handle complicated Hamil- 
tonians is through Matrix Product Operators (MPO). 
MPO's provide a convenient way to reason about op- 
erators in MPS algorithms and have been employed ex- 
tensively in time-dependent MPS simulations. [31] Here 
we provide a brief analysis of this approach for quantum 
chemistry. The basic idea in an MPO is extend the ma- 
trix product like representation to operators. Writing the 
Hamiltonian in the occupation number basis as a general 
rank-2A: tensor, [^]^i- -^/e^i - "^fc^ it is decomposed into a 
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FIG. 4. Graphical representations of MPO. The Hamiltonian 
operator is a rank-2A; tensor (top panel). This can be divided 
into a contracted product of site-independent tensors similarly 
to an MPS, leading to a Matrix Product Operator represen- 
tation (MPO, bottom panel), in which each site tensor is a 
rank-4 tensor. 

set of tensors analogously to an MPS, 

This is illustrated graphically in Figure 2] where W^*] is 
a rank-4 tensor and the contraction of horizontal edges 
describes the quantum "entanglement" of the Hamilto- 
nian operator. To decompose H exactly, the dimension 
of the horizontal edge of W^^^ needs to be exactly the 
same as the number of complementary operators, that 
is 0{k'^). Consequently, the cost of computing the ex- 
pectation value with an MPO representation 
becomes 0{M^k^ + M'^k^). Note that this cost is larger 
than in the complementary operator approach. The dif- 
ference arises because we have not considered the sparsity 
of the individual W^^^ tensors that arise in the Hamil- 
tonian decomposition. However, incorporating element- 
wise sparsity into an MPO algorithm eliminates much of 
the conceptual and algorithmic simplicity of the MPO 
approach. Consequently, in our view, the complemen- 
tary operator algorithm is a more practical and efficient 
route for MPS computations with quantum chemistry 
Hamiltonians. The relative benefit of using complemen- 
tary operators versus tensor product operators (TPO's) 
is even greater for TTNS than for MPS, thus we focus on 
the complementary operator approach when considering 
trees. 



TREE TENSOR NETWORK STATES (TTNS) 

Tensor Network States (TNS) are mathematical gen- 
eralizations of the MPS that can code more general en- 
tanglement networks. The form of a TNS wavefunction 
is directly analogous to an MPS wavefunction, 

|*)= J2 ttr[A"i -...A"- ■...A"'=]|ni...nfe), (15) 

ni...nfc 




FIG. 5. Graphical representation of TTNS. Left panel shows 
overall structure of TTNS as described in Eq. ([T5|), and right 
panel shows the one-site wavefunction spanned by the renor- 
malized basis, as described in Eq. (|16p . Note that physical 
indices (vertical bonds) are omitted in the right panel. 

the only difference being that A"^* is now a tensor, rather 
than a rotation matrix as in the MPS, and the multipli- 
cation operator • together with ttr denotes a general con- 
traction over tensor indices. The flexibility of TNS wave- 
functions arises from the fact that whereas the matrices 
in an MPS can only be contracted along a ID-lattice, 
there are many different ways to connect general tensors 
together to form a network of entanglement. 

Tree tensor network states (TTNS) are a special class 
of TNS where the tensors are connected as a tree as 
shown in Figure [5l A tree is a graph that has no loops, 
which leads to many simplifying mathematical properties 
that parallel those of a MPS (indeed, an MPS is simply 
a tree with Z = 2 legs). For example, at a given site i in 
the tree, we can define renormalized bases for each of the 
Z legs connected to the site. This allows use to rewrite 
(fT5l) as 

i^)= E |6^..6f n.) (16) 

bl...bfn^ 

where \bf) is the renormalized basis in the a-th branch 
of site i. This basis is defined by recursively contracting 
tensors in the branch from the leaves up to site i. Thus 
\bf) is obtained as 

1^?)= E <l..,f-^,?\'^-'^f~'^^) (17) 

3 3 -J 

where sites j are adjacent to i in the branch. Note that 
for a Z degree TTNS with M renormalized states, the 
tensor A^^ b^'^b"^ 0{M^) elements. 

3 "' 3 i 

Analogous to the rotation matrices in MPS, the tensors 
A^^ in a given branch around site i can be chosen to 
satisfy orthonormality constraints, rendering the TTNS 
in canonical form, 

b^....b^~^ni 

3 3 ^ 

As a result, the renormalized basis states \bf) are or- 
thonormal. 
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The above mathematical properties make a DMRG en- 
ergy optimization algorithm for TTNS very similar to 
that for MPS. Similarly to the DMRG algorithm for 
MPS, we optimize one site at a time. The TTNS is ex- 
pressed in canonical form around site z, then the coeffi- 
cient tensor -0^^ is optimized, and the sites of the tree 
are traversed during the sweep. The computational chal- 
lenge is once again how to efficiently compute the repre- 
sentation of the Hamiltonian in the renormalized bases, 
namely {bj ...bf ni\H\b'] ...b'f n' i) ^ and its action on the 
coefficient vector 1/?"^% Hi/j. As we discussed for the case 
of MPS, for quantum chemistry Hamiltonians the com- 



plementary operator approach is most natural. 

We rewrite the Hamiltonian in complementary opera- 
tor form by partitioning into Z -\- 1 blocks. A, B, C, D, 
and so on, corresponding to the Z branches around site i, 
and site i itself. The Hamiltonian is re-expressed in terms 
of operators acting on each of the blocks separately, 

^ ^ ^ 

+ Habc + Hacd + Hbcd + ••• 

+ Habcd + ••• 



where, 

I 

Hab= E (af^f + af^/ + afV^^^ (20) 
ieAjeB i>jeA 



habc= E + «f 4^0 + ^ («^?^^f + «f^«?4^0 

ieAjeB ieAjec 



ieBjec 



Habcd = ^ ^ VijikCtf^af^a^aP + permutation. (22) 
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ieAjeB 
kec,ieD 



As in the case of MPS, the full matrix representation 
of H is never built explicitly as the storage requirements 
would be immense. Instead only Hip is computed in the 
Davidson algorithm. Note that when computing Htp the 
order of multiplication of the various terms in the com- 



plementary operator decomposition of H is important. 
For example, the term afd^P!^^ip in which A and C are 
large blocks with 0{k) sites and M renormalized states 
(indices bi and 62 below) and B is the site being opti- 
mized, with 4 states (index rii below) should be computed 
as 



ieAjeB 



[af] 



ieA,bi \b2,ni 



KjeB 



2 02 



(23) 



In the above form, T^jeBi^f] ^ [^j^] is of 0{WM^k) 
cost, is of (9(64M^+ifc) cost, and finally 

is of (9(4M^+i/c) cost. 

We have carefully optimized the order of tensor con- 
tractions for each of the terms in (|2Q|) - ([22]) . obtaining the 
complexities shown in Table HI For a general Z degree 
tree, the total cost of a single Htjj (sigma vector) compu- 
tation is 0{M^^^k^ + M^k^) per sweep. To this must 
be added the construction of the matrix representations 



of the components of H in the different blocks (renor- 
malization steps), which costs 0{M^^^k^ + M^k^) per 
sweep (The renormalization is described explicitly in Ap- 
pendix A). In the case of Z = 2 (MPS) we have already 
noted that the cost per sweep is lower than in the general 
case {0{M^k^ + M^k^)). This is due to the absence of 
several terms in Table HI The same is true for the Z = 3 
tree, where the total cost of a sweep is 0{M^k^ + M'^k^) 
(Hi/j) and (9(M^/c^ + M^/c^) (renormalization). Because 
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TABLE I. Complexity of optimal tensor contractions for H^ip 
per site. Multiplying by 0{k) gives complexity per sweep. 



(A) 



terms in H 

Qj^ Qjj Qij 

a- 'a. 'a^ai 



complexity 




I 4M k ^rwT>r 

left block 1 i+V 



FIG. 6. Half-Renormalization (HR) algorithm on TTNS. (A) 
One-site algorithm: Z — 1 system blocks {|&?)} are mapped 
into one system block. (B) Two-site algorithm: additionally, 
Z — 1 environment blocks {|^f+i)} are mapped into one en- 
vironment block. The half renormalized block contains only 
AM states for any Z > 2. 



of the special efficiency of the Z = 3 tree, our later com- 
putations focus on this kind of tree. 



Half-renormalization 

We have found additionally that it is possible to sig- 
nificantly reduce the computational prefactor of a TTNS 
DMRG calculation through an additional step we call 
half-renormalization (Fig. [6j). Half-renormalization in- 
volves first constructing an exact mapping of the TTNS 
onto an equivalent Z = 2 MPS, then carrying out the 
Hip operations in this simpler representation. Although 
the mapping is itself expensive (and retains the full com- 
putational scaling of the TTNS sweep described above), 
it need only be carried out once per site, while the Hip 
operations typically need to be carried out many times 
per site during a Davidson diagonalization. To map a 
TTNS onto an MPS, we consider the coefficient tensor 
at site i, ^z- Through an SVD, this tensor can be 




(B) 



FIG. 7. Sweep algorithm on a tree (Z = 3), by depth- first 
search with backtracking, where the labels indicate the order 
of the search. (A) One-site algorithm: starting from the cen- 
ter site, 1 sweep contains 18 one-site optimization steps, and 
(B) Two-site algorithm: starting from the center site and one 
adjacent site, 1 sweep contains 12 two-site optimization steps. 



exactly decomposed into a rank-3 tensor and a residual 
tensor, 

.ni(TTNS) _ V^rr o T/n^f 



i(MPS) 



(24) 



where U is the residual tensor and SV^ is an MPS-like 
pseudo-coefficient tensor tp^^2^^^\ Note that 6* is a di- 
agonal matrix with only 4 x M non-zero singular values, 
thus ip^f^z^^^^ has O(M^) values, similar to an MPS co- 
efficient tensor with M renormalized states. /J^i . ^z-i^ 
defines the half-renormalization mapping from the states 
on Z—1 branches, a = 1 . . . Z—1 of a general tree to a sin- 
gle set of 4M renormalized states \bf) on a single effective 
branch. These new "half-renormalized" states, together 
with the states on the remaining branch a = Z, define 
the left and right states of an effective MPS, together 
with the pseudo-coefficient tensor tp^^2^^^\ In the case 
of Z = 3, the mapping costs 0{M^k^ -\- M'^k^) per sweep 
which leads to the same scaling as the one-site algorithm 
without half-renormalization. However, the HiIj oper- 
ations now carry the MPS cost, namely 0{M^k^) per 
sweep, leading to a significant savings in computation 
time. As the mapping is exact, there is no approxima- 
tion involved, although one drawback is that we typically 
see slower convergence, as at a given site we are optimize 
only the pseudo-coefficient tensor rather than the full co- 
efficient tensor of the TTNS. 

The full DMRG sweep on a tree is carried out using 
depth-first search with backtracking. Applying this to 
a Z = 2 tree, as in an MPS, recovers the usual DMRG 
sweep algorithm. Figure [71^ A) shows the details of the 
one-site sweep algorithm on TTNS. 

(1): Define any site as the root (depth = 0). Construct 
the TTNS in canonical form at the site, construct 
the renormalized states and operators by contract- 
ing from the leaves to the root site. 
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FIG. 8. Block structure of the two-site algorithm. Top panel 
shows the two-site MPS which contains one system block, 
system site, environment site, and one environment block. 
Bottom panel shows the two-site TTNS which contains Z — 
1 system blocks, system site, environment site, and Z — 1 
environment blocks. 

(2) : Optimize the wavefunction at site i with David- 

son diagonalization. Renormalize to construct the 
canonical form at a neighbouring site. 

(3) : Continue (2) by carrying out a depth- first search 

with backtracking to optimize the entire tree. 

(4) : Continue (2)- (3) until the energy is converged to a 

target accuracy. 

Two-site TTNS algorithm 

Although our discussion of MPS focused for simplic- 
ity on the one-site MPS and the corresponding one-site 
DMRG algorithm, it is well known from practical experi- 
ence that one-site DMRG calculations tend to get stuck 
in local minima and suffer from poor convergence charac- 
teristics. Thus, the two-site MPS and DMRG algorithm 
are more commonly used.0, [lH, [36| In the case of MPS, 
the two-site MPS wavefunction is obtained by modifying 
the coefficient tensor to span two sites, 

l*)= E Cmt J^i-i^i^i+i'-i+i) • (25) 

A Z-branch two-site TTNS is modified in a similar way 
E C'^^-.^ ....-^ |6,^...6f-in,n,+i6l+,...6f+li) 

nirii+i * * 
hihi+i 

(26) 

An important difference however, between a general two- 
site tree and two-site MPS is that whereas around sites 
i, i + 1 in an MPS we can define left and right renormal- 
ized basis, just as for an one-site MPS, in the case of a 
TTNS, there are 2Z — 2 branches around sites + 1 
(Fig. [8j), rather than the Z branches around a single- 
site. Consequently, the naive cost of the Htp opera- 
tions in a two-site TTNS DMRG sweep using complemen- 
tary operators is much higher than that for an one-site 
TTNS sweep, with a prohibitive cost of ©(M^^-^A:^ + 



^2Z-2^5j However, by employing two sets of half- 
renormalization steps, we can map the two-site TTNS 
onto a two-site MPS, reducing the cost of the Htp opera- 
tions to only 0{M^k^) per sweep with an additional half- 
renormalization cost of 0{M^^^k^ + M^k^) per sweep 
for general Z degree trees, or 0{M^k^ ^M'^k^) per sweep 
in the case of Z = 3. This is the same cost as an one- 
site TTNS sweep with half-renormalization, and thus the 
two-site TTNS sweep becomes practical. We compare the 
cost and convergence characteristics of the one-site and 
two-site TTNS sweeps with half-renormalization in our 
later calculations. 

Tree Shape and Site Ordering 
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FIG. 9. Approximate degree- fixed minimum spanning tree 
(MST, Z = 3) for 24 orbitals from a RHF/cc-pVDZ calcu- 
lation of the water molecule. Left panel shows a representa- 
tion of the exchange integral Kij^ where solid lines denote 
Kij > 0.10, dashed lines denote Kij > 0.07, and dotted 
lines denote Kij < 0.07. Colored lines are connected lines 
in the right panel and gray lines are ignored interactions. For 
Kij < 0.07, only connected lines are shown. Labels in the 
right panel indicate the MOs (as indexed by energy). 

When carrying out a quantum chemistry DMRG calcu- 
lation using an MPS, it is necessary to choose a mapping 
of the sites to the ID lattice. Generally, this should be 
done to minimize entanglement between distant sites on 
the lattice, js^-Q but computing the entanglement and 
carrying out an exact minimization are costly procedures. 
In practice, an approximate proxy for the entanglement 
between orbitals is constructed, and it is approximately 
minimized. j3] One example of such a proxy which has 
been used in prior DMRG studies is a weighted exchange 
integral KijDfj^ where K^j is the exchange integral be- 
tween orbitals i and j, Dij is the separation on the lattice, 
and n is an adjustable parameter. 

In the case of ordering for trees, we need to consider 
not only the mapping of the orbitals onto a given tree, 
but also the shape of the tree as well, even if we restrict 
ourselves to trees of a fixed maximum degree. Since the 
computational cost of a sweep also depends on the tree 
shape (changing the prefactor associated with the cost of 
the sweep) the tree which gives the lowest energy for a 
given M is not necessarily the most efficient tree to use 
in practice. 
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TABLE II. The ground state energies of water molecule calculated with two different trees. 



M 




MST 




MET 


E 1 a.u. 


CPU time / sec. sweep 


E 1 a.u. 


CPU time / sec. sweep 


100 


-76.242805 


307.9 


-76.242242 


220.6 


200 


-76.243652 


1209.9 


-76.243491 


826.5 




FIG. 10. Minimum-entangled tree (MET, Z = 3) for 
24 orbitals from a RHF/cc-pVDZ calculation of the water 
molecule. Left panel shows the stepwise construction of MET. 
First, 24 sites are divided into 1 center site and 3 branches 
containing 7, 8, and 8 sites. Second, 7 sites are divided into 1 
site and 2 branches containing 3 sites for each, and 8 sites are 
divided into 2 sites and 2 branches containing 3 sites for each. 
Finally, we used a genetic algorithm to map the orbitals to 
the tree sites as shown in right panel. Labels in right panel 
indicate MO indices (as indexed by energy). 



To illustrate these issues we consider the water 
molecule in a cc-pVDZ basis set, which has 24 orbitals. 
We have considered two shapes of trees: a degree-fixed 
Minimum- Spanning Tree (MST) shown in Figure [H and 
a Minimum-Entangled Tree (MET). The MET is defined 
as the shape of tree, for a given degree, where the num- 
ber of renormalized states required to achieve an exact 
calculation is minimized. Its construction is shown in 
Figure [TOl Because of its balanced nature, the MET also 
minimizes the prefactor of the cost of the TTNS calcula- 
tion for a given M. 

For each of these trees, we mapped the orbitals onto the 
sites by a genetic algorithm [40] that minimized the cost 
function KijD^j^ where Dij is the counting distance 
between the sites i and j in the tree. 

The ground state energies computed with the two dif- 
ferent trees and two different M (100 and 200) are sum- 
marized in Table ITTl From these calculations, we see that 
the MST gives a slightly better energy than that the 
MET, but the MST sweep is 50% slower than the MET 
sweep for the same M, due to the unbalanced nature of 
the MST. For larger M, because the MET minimizes the 
number of renormalized states for the exact calculation, 
we expect it to eventually give a lower energy for a given 
M than the corresponding MST. Although these results 




1 sites (g = 2) 22 sites (g = 3) 46 sites (g = 4) 



FIG. 11. Hydrogen atoms on Cay ley-trees (Z = 3, ^ = 2(10 
sites), 3 (22 sites), and 4 (46 sites)). Top panel shows the 
actual structures of the hydrogen trees, in which red atoms 
are the core region, orange atoms are for g — 2, yellow atoms 
are for ^ = 3, and green atoms are for g = 4. The distance 
between adjacent hydrogens is 2.0 Bohr and torsional angle 
between different generations is 30 degrees. Bottom panel 
shows corresponding Cayley-tree diagrams. The red-dotted 
line denotes the site ordering used in the MPS. 

are system dependent, they indicate the importance of 
tree shape in determining the cost of practical TTNS cal- 
culations. We adopt the MET in our calculations unless 
another tree shape is trivially indicated by the molecu- 
lar structure, e.g. in a tree shaped molecule such as a 
dendrimer. 



ILLUSTRATIVE CALCULATIONS 

We now present several benchmark applications to 
molecular systems to understand the performance of 
TTNS and to compare with that of MPS. For the MPS 
calculations, we used our TTNS code with Z = 2. This 
was to allow a fair comparison of timings using the same 
implementation. 

To start, we consider an idealized system: hydrogen 
atoms on Cayley-trees in an (orthogonalized) minimal 
ST0-3G basis. This is a model system where TTNS are 
expected to work very well. To avoid nearly overlapping 
hydrogens which would arise in a planar geometry, the 
tree structures are taken to be slightly twisted as shown 
in Figure [TTJ All geometries are provided in the supple- 
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^ 1 .OE-03 - 



HR 2-site (converged at 16 cycles) 
21 1.3 sec. 



Full 2-site (converged at 4 cycles) 
3468.2 sec. 

1 1 

5 10 

sweep iterations 



TABLE III. CPU times per sweep in sec. for MPS and TTNS 
calculations for the nitrogen dimer using (lOe, 26o) active 
space and cc-pVDZ basis sets. 



FIG. 12. Energy convergence and total CPU time com- 
puted with normal two-site TTNS (denoted Full) and half- 
renormalized two-site TTNS (denoted HR-TTNS). 



mentary information. [45| 

We first illustrate the importance of the half- 
renormalization algorithm by comparing the normal two- 
site TTNS algorithm and the HR-TTNS two-site algo- 
rithm on a small 10 site hydrogen tree. The energy con- 
vergence along the sweeps and total CPU time are shown 
in Figure [121 Though the convergence per sweep is slower 
in the HR-TTNS algorithm, the total CPU time is much 
smaller than in the normal two-site TTNS algorithm due 
to the much less expensive Hip operation. Thus, the re- 
maining calculations have been done using the HR-TTNS 
algorithm, which we henceforth shorten to TTNS. 

Next, we compare the performance of the MPS and 
TTNS on three hydrogen trees of different sizes: 10 sites, 
22 sites, and 46 sites. We focus on two aspects: energy 
convergence with M, and CPU time per sweep. In gen- 
eral, the energy convergence of the TTNS is much faster 
than the MPS, both with increasing M and both as a 
function of CPU time, as seen in Figure [13] for the 22 site 
tree. (In the 10 site tree, the strange behaviour of the 
CPU time, where larger M required smaller CPU time, 
is an artifact of our Davidson diagonalization implemen- 
tation, which required a large number of iterations when 
M was very small). In the 46 site tree, we only carried 
out calculations with small M because of the large com- 
putational cost. For the CPU time per sweep, we found 
that a calculation with MPS with M = 200 (6182 sec. / 
sweep), and TTNS with M = 100 (6839 sec. / sweep) 
were comparable. However, the corresponding energies 
were computed to be -22.705014^^ and -22. 9582 19^;/,, 
respectively and thus the TTNS gave a much better en- 
ergy than the MPS at a comparable computational cost. 
These results demonstrate that TTNS are more cost ef- 
fective than MPS when the system is tree-shaped. More- 
over, the relative benefits of the TTNS increase as the 
tree size increases. 

We now consider a model system designed to mimic 
more difficult molecular structures for MPS and TTNS. 
In the previous tree structures we used a local atomic ba- 



R/ A 


MPS 


TTNS 


M = 500 M = 1000 


M = 250 M = 500 


1.1208 


706 2998 


794 5088 


1.4288 


1027 4684 


1337 7510 


1.9050 


1155 5053 


992 5945 



sis. This, together with the underlying tree connectivity 
of the molecule, allowed the TTNS to completely exploit 
the local nature of correlation in the system. Though 
the system was not linear, the MPS still benefitted from 
the local basis in the calculations, as not every orbital 
was correlated with every other. (For similar reasons, in 
DMRG calculations localized molecular orbitals are often 
used to minimize long-range entanglement). However, in 
some situations, it may be unavoidable for an MPS or 
TTNS to describe some long-range entanglement, either 
due to the underlying physics, or due to a poor mapping 
of the molecular structure onto a chain or tree. To mimic 
this situation we carry out MPS and TTNS calculations 
in a canonical molecular orbital basis. Because these or- 
bitals are delocalized, every orbital can be considered to 
interact with every other. We choose as our model sys- 
tem a set of hydrogen chains in a minimal ST0-3G basis. 
Although these are chain molecules, the use of canonical 
molecular orbitals means that this is not an ideal system 
for the MPS. 

Figure [14] shows the energy convergence as a function 
of M and CPU time per sweep. Both MPS and TTNS 
exhibit very slow energy convergence as a function of M 
due to the use of canonical molecular orbitals, but it is 
clear that the convergence with M is always better in 
TTNS than in MPS by a factor of 2 or more. As a func- 
tion of CPU time, we find that the decreased M means 
that TTNS performs better in comparison with MPS, al- 
though the improvement is slight due to the higher scal- 
ing of TTNS with M. 

These results provide some promise that in general 
molecules without a tree or chain structure, TTNS can 
perform better than MPS, although this involves a deli- 
cate balance between the decreased number of renormal- 
ized states used by the TTNS, and the higher computa- 
tional scaling with M. To examine this in a more realistic 
setting, we now consider calculations on two benchmark 
diatomics, the nitrogen dimer and the chromium dimer. 

The bond dissociation curve of nitrogen dimer is often 
used as a good benchmark to evaluate whether a method 
can describe strong electron correlation correctly. We 
evaluated ground state energies at three points on the 
bond dissociation of nitrogen dimer, 1.1208(i?e), 1.4288, 
and 1.9050 A. A frozen core active space (lOe, 26) with 
a cc-pVDZ basis set was employed as used in previous 
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FIG. 13. Energy convergence of hydrogen trees plotted with respect to the number of renormahzed states M and CPU time 
(sec.) per one sweep; (A) and (B) show energy versus M and CPU time per sweep, respectively, for 10 sites, and (C) and 
(D) show energy versus M and CPU time per sweep, respectively, for 22 sites. All calculations are for the triplet ground state 
energy in an orthogonalized ST0-3G basis. 



DMRG calculations and full-CI calculations [4lj (see 
supplementary information for the tree graphs and the 
site orderings|45|). Figure [15] shows the energy errors 
from full-CI results [41^ as a function of bond length and 
CPU times per sweep are summarized in Table lllli In 
TTNS, half the M can be used at = 1.1208 and 1.4288 
A and a quarter the M can be used at = 1.9050 A as 
compared with MPS. TTNS with half the M of the cor- 
responding MPS required approximately the same CPU 
time per sweep. This indicates that the TTNS gives 
comparable and/or slightly better performance at short 
bond-lengths but much better performance at long bond- 
lengths compared to the MPS in this molecule. We con- 
clude that the TTNS works better in the case of the ni- 
trogen dimer than the MPS. 

Since the chromium dimer has an unusual multiple 
bond, its ground state is very complicated and difficult 
to describe by conventional methods. Recently, DMRG 
calculations of the chromium dimer have been carried 
out for relatively large active spaces. [9|,|42| We performed 
MPS and TTNS calculation at 1.5 A usins: the same ba- 
sis sets and the same active space (24e, 30o) as in earlier 
DMRG calculations [42] (see supplementary information 
for the tree graph and the site ordering [45j). Figure [16] 
shows the energy convergence as a function of M and 
CPU time per sweep. Since the previous DMRG calcu- 



lation employed very large M (up to 10000), the energy 
reported here is far from convergence because we only 
employed M up to 150 in the TTNS. Nonetheless com- 
pared to our MPS benchmarks, we find that smaller M 
can be used in the TTNS. In comparing CPU time per 
sweep, however, we see that MPS gives much better per- 
formance than the TTNS in this molecule. Thus, we 
conclude that the MPS works better in the case of the 
chromium dimer than the TTNS. 

These two benchmark calculations on the nitrogen 
dimer and chromium dimer show that the performance 
of the MPS and TTNS in general molecules depends sen- 
sitively on the electronic structure and the nature of the 
quantum entanglement of the molecule, thus their rela- 
tive merits must be determined on a molecule by molecule 
basis. 

Finally, to demonstrate the power of TTNS, we con- 
sider a TTNS calculation on more realistic tree-shaped 
molecules. Stilbenoid dendrimers are prototypical tt- 
conjugated dendrimers, with each unit being a stilbene 
fragment. These dendrimers are potentially attractive 
for chemistry because photo-induced electron-transfer or 
exciton-transfer may proceed from the leaves to the core, 
mimicking a biological antenna system (we note that den- 
drimer systems have also been previously studied using 
semi-empirical DMRG [Isl, Q). Although it is the ex- 
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FIG. 14. Energy convergence of hydrogen chains with CMOs plotted by number of renormahzed states M and CPU time (sec.) 
per one sweep; (A) and (B) showed those of M and CPU time per sweep, respectively, for 20 sites, and (C) and (D) showed 
those of M and CPU time per sweep, respectively, for 30 sites. Singlet ground state energy was computed with ST0-3G basis 
sets. 
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FIG. 15. Energy errors for bond dissociation of nitrogen dimer 
using (lOe, 26o) active space and cc-pVDZ basis sets. Refer- 
ence energies were from previous full- CI work in Ref. [411 . 



cited states and dynamical properties of these systems 
that are of primary interest, here we focus on the ground- 
state energy for benchmarking and reserve the study of 
excited states to future work. 

We consider three different sizes of stilbenoid den- 
drimers denoted by the generation g as shown in Fig- 
ure [m We carried out 7r-full valence MPS and TTNS 
calculations with ST0-3G and 6-31G basis sets. Single 



valence calculations were performed for g = 0, 1, and 
2 with (14e, 14o), (46e, 46o), and (llOe, llOo) active- 
spaces, respectively, and double valence calculations were 
performed for ^ = and 1 with (14e, 28o) and (46e, 92o) 
active-spaces respectively. To construct the tree graph, 
TT-orbitals computed from a RHF calculation were local- 
ized for occupied and unoccupied spaces separately. Lo- 
calized MOs were grouped for each ethylene and ben- 
zene fragment, and were ordered on the tree according 
to the underlying dendritic structure, as shown in Fig- 
ure [iHl The orderings within each fragment were de- 
termined to put strongly interacting pairs (evaluated by 
Kij) on neighboring sites. 

For the ^ = dendrimer (stilbene), the MPS and 
TTNS gave similar energy convergence as a function of M 
in the single- valence (ST0-3G) calculation, and the MPS 
gave somewhat better performance with respect to CPU 
time per sweep (see Figure [T9R andfTOB). In the double- 
valence (6-3 IG) calculation, the TTNS gave an improved 
energy convergence with M compared to MPS, and the 
MPS and TTNS gave almost the same performance with 
respect to CPU time per sweep (see FigurefTOD andflOD). 
The competitive performance of MPS relative to TTNS 
in the ^ = dendrimer reflects the very small size of the 
system, which is almost linear in nature and thus nearly 
ideal for MPS. 
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1. OE-01 




FIG. 16. Energy convergence of the chromium dimer at 1.5 A using a (24e, 30o) active space and Ahhichs' SVP basis set. (A) 
Plotted as a function of M and (B) as a function of CPU time per sweep. Reference energies were from previous DMRG work 
in Ref . Si. 




g=0 g=1 g=2 

FIG. 17. Structures of stilbenoid dendrimers ^ = 0, 1, and 2. 
Geometries were optimized at the B3LYP/cc-pVDZ level of 
theory. 



For the larger g = 1 dendrimer, in the single- valence 
space the TTNS required a quarter the M to obtain the 
same energy accuracy as the MPS. The computational 
cost of the TTNS for this accuracy was also much lower 
than that of the MPS (see Figure HHP and HHP). The 
double- valence calculations were very expensive, hence 
we only performed calculations with M = 100 for MPS 
and with M = 50 for TTNS. Although the resulting en- 
ergies were not converged to chemical accuracy, the cor- 
relation energies {E - E^f) of the MPS, -0M4736Eh, 
and the TTNS, -0.67S678Eh, indicate that the TTNS is 
once again much more accurate. 

Finally, for the g = 2 dendrimer, even the single- 
valence active space calculation was quite expensive 
with our implementation, consequently, we only per- 
formed single- valence calculations using M = 100 for 
MPS and M = 50 for TTNS. These calculations had 
comparable timings. The computed correlation ener- 
gies were -1.943651£;^ and -2.021378^;/, with MPS and 
TTNS, respectively, indicating once again that the TTNS 
works much better than the MPS in this large den- 
dritic molecule. Although we could not reach fully con- 
verged energies in this work, the ability to even approxi- 
mately target such systems with TTNS demonstrates the 



promise of the technique for complex systems. 



CONCLUSIONS 

In this work we investigated tree tensor network states 
(TTNS) for quantum chemistry. We formulated an ef- 
ficient tree tensor network algorithm that is analogous 
to the density matrix renormalization group (DMRG) 
algorithm in quantum chemistry for matrix product 
states (MPS). We introduced the additional step of 
half-renormalization that greatly reduced the computa- 
tion cost. We found that our TTNS calculations were 
competitive with MPS and DMRG calculations in gen- 
eral molecules, requiring significantly fewer renormalized 
states for the same accuracy, although this did not al- 
ways translate into a savings in computational time. In 
tree like molecules, TTNS were clearly superior to MPS 
requiring both fewer renormalized states and less time 
to reach the same accuracy. This bodes well for the ap- 
plication of TTNS to study a wide class of interesting 
optically active systems based on dendritic structures, as 
illustrated in our calculations on stilbenoid dendrimers. 

Tree tensor network states are one of the simplest gen- 
eralizations of the matrix product states, because the 
tree network structure has no cycles. More complex ten- 
sor networks which describe even more general entangle- 
ment require the treatment of cycles. The improvements 
observed with tree tensor networks here suggests that 
investigating these more complex classes will be fruit- 
ful, particularly to describe quantum chemistry of larger 
systems, where the molecules or materials acquire an ex- 
tended two-dimensional and three-dimensional structure. 
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ethylene frag. 


benzene frag. 




FIG. 18. Tree graphs for TTNS calculation of stilbenoid dendrimer = 1). (A) Selected molecular orbitals of the localized 
ethylene and benzene fragments, (B) Tree graph of the single-valence (ST0-3G) calculation, and (C) tree graph of the double- 
valence (6-31G) calculation where the filled circles represent local tt and tt* orbitals and empty circles represent local 3p orbitals. 
Note that identical orbitals are omitted in the tree graphs. 
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(28) 



APPENDIX A: OPTIMAL TENSOR 
CONTRACTION FOR TTNS 
RENORMALIZATION 



During the TTNS renormalization procedure, it is nec- 
essary to construct matrix representations of operators 



act on b}...b, 



z-i 



rij 



Hi) \bf~^) ... \bl). Because of the 



quartic terms in the Hamiltonian, this can involve prod- 
ucts of up to four operators. The corresponding product 
formulae follow. 

The matrix representation of an operator Oa which 
acts only on renormalized states in block a is given by. 



]^Q^l^b']...b'f-^n'ibl...bf-^ni 

= {b^]\...{b^f-'\{n\\d^\n,)\bf-')...\bl) 



(27) 



[da 
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xP{{b't^'\...{bT'\{n'^\.Oa) 



xP{{bT'\-{b't'\{n'i\,dc,) 
xP{{b'^^'\...{b'f-'\{n\\,60) 
xPi{bT'\-{b'f-'\{n'.\,d,), 



[ds 



xP((C+V..(&f"'lR|,6«) 



xP{{b'^^'\...{b'f-'\{n'i\,d0) 
xP{{b'r'\-{b'f-'\{n'i\,d,) 



xP{{bT\-{b'r \{n'^\,Os). 



(29) 



(30) 



where [Oa]^ * is the matrix representation of Oa 
spanned by states \bf) and P is a parity operator which 
gives +1 or —1 depending on the particle numbers of 
{b'f~^^ ...b^f~^n'i\ and Oa- Similarly, those of operators 
Oq;0/3, OaOjsO^^ and OaOjsO^Os are computed as fol- 



Here, the representation of operator 

[0]6'1...5'f-^n%5l...6f-^n, ^ /^j^Z ^ ^j^Z j^atrix. 

In the renormalization step, this is reduced to M x M 
matrix [O]^'^ , spanned by renormalized states \bf) 
using a tensor A"^^ as. 
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FIG. 19. Energy convergence of stilbenoid dendrimers ^ = and g — 1 plotted by number of renormalized states M and CPU 
time (sec.) per sweep; (A) and (B) shows energy versus M and CPU time per sweep, respectively, for ^ = 0, single valence 
(14e, 14o), (C) and (D) shows energy versus M and CPU time per sweep, respectively, for ^ = 0, double valence (14e, 28o), 
(E) and (F) shows energy versus M and CPU time per sweep, respectively for g — 1^ single valence (46e, 46o). All energies are 
for the singlet ground-state energy. 
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(31) 

It is necessary to minimize the cost of tensor contrac- 
tions for TTNS renormalization (|3T]) is just as for Hip^ 
as described in the text. To compute representations of 
complementary operators, each complementary operator 
can be further divided into Z — 1 renormalized blocks and 
the site i itself. We carefully minimized these tensor con- 
traction costs for each complementary operator and their 



complexities per site are summarized in Table IIVI for a 
general tree. Thus, for a general tree, the computational 
complexity of a TTNS sweep is 0{M^^^k^ + M^k^). 
As discussed in the text, the complexity can be further 
reduced to 0{M^k^ + M'^k^) for Z = 2 (MPS) and 
0{M^k^ + W^k^) for Z = 3. 
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